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Temperature and density extrapolations in canonical ensemble Monte Carlo 

simulations 

A. L. Ferreira and M. A. Barroso 
Universidade de Aveiro, Departmento de Fisica, 3810-193 Aveiro, Portugal 

We show how to use the multiple histogram method to combine canonical ensemble Monte Carlo 
simulations made at different temperatures and densities. The method can be applied to study sys- 
tems of particles with arbitrary interaction potential and to compute the thermodynamic properties 
over a range of temperatures and densities. The calculation of the Helmholtz free energy relative to 
some thermodynamic reference state enables us to study phase coexistence properties. We test the 
method on the Lennard- Jones fluids for which many results are available. 

PACS numbers: 05.10.Ln, 05.20JJ, 64.70.Fx 
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I. INTRODUCTION 

Histogram and multiple-histogram methods have been 
proposed as an optimized way of analyzing Monte Carlo 
data pi pf. These methods can be included in the more 
general class of reweighting methods g . The idea is to 
combine a given set of standard Monte Carlo simulations 
to get improved estimates of observables in a given pa- 
rameter region. A related idea is to sample a suitably 
chosen probability distribution rather than a given sta- 
tistical mechanics ensemble. The sampling distribution is 
such that the configuration space visited is typical of the 
interval of thermodynanric paranreters of interest thus 
allowing the reconstruction of the appropriate statistical 
mechanics ensemble. The methods of umbrella sampling 
[grff, multicanonical [|]|| and expanded ensemble meth- 
ods [J15],[L6| can be seen as belonging to this class. 

In this article we show how to combine NVT Monte 
Carlo simulations made at different temperatures and 
volumes. Our method is a generalization to volume ex- 
trapolations of the multiple histogram method. We fur- 
ther show that the method can be applied not only to sys- 
tems of particles that interact through interaction poten- 
tials that have a simple scaling with particle distance but 
also to those with arbitrary distance dependence. Fur- 
thermore, as we are able to calculate relative free energies 
as a function of volume and temperature, the method can 
be applied to study phase coexistence properties [gJQ]- 

Several simulation methods have been proposed to 
study phase coexistence properties. In the Gibbs en- 
semble Monte Carlo |1(J two simulation boxes equili- 
brate by exchanging particles and volume and the sys- 
tem separates into two phases, each one located in one of 
the boxes. Grand-Canonical ensemble simulations with 
multiple histogramming have been used to study critical 
properties and finite-size scaling in fluid systems p2|-[l4[ . 
However, due to the low probability of particle exchanges 
or insertions at high densities these methods cannot be 
used to study dense phases. For such systems special 
methods have been suggested that rely on the calcula- 
tion of the absolute free energies of the two phases |L7J . 
Recently a new method based on Gibbs-Duhem integra- 



tion was proposed and the concept of pseudo-ensenrbles 
was introduced |19-21|. For all of these methods the use 
of the multiple histogram technique can be a valuable 
auxiliary tool [ p7pl| ] . Our method can also be applied to 
the study of solid fluid coexistence p8[ . 



II. THE METHOD 

Consider a system of N interacting particles contained 
in a box of volume Vq. For sinrplicity we consider a pair- 
wise additive interaction potential and {f^} denotes a a 
given configuration of particle coordinates. The total po- 
tential energy of the system in a configuration is given by 
E({?i}) — J2<i j> u {\ri~- rj\), where the sum runs over all 
pairs of particles. Uniformly expanding the system from 
the volume Vq to the volume V changes the configuration 
from {fi} to {fit}, such that f,/ = (V/Vq) 1 ^ 3 ^. The en- 
ergy of the system of volume V in the new configuration 
is given by, £({?«/}) = E< lJ > «((^/Vo) 1/3 |fi - r* |). 

We will show that it is always possible to find a set of n c 
variables, C„({ri}), with < n < n c — 1, that depend on 
the particle coordinates. These variables can be seen as 
coordinates of a column vector, C = (Co, C±, ..., C„ e _i). 
Their choice is arbitrary provided two properties are ful- 
filled. First, it should be possible to write the potential 
energy in terms of these variables. Second, there is a 
known linear relation between the value of the variables 
in the expanded system of volume V and their value for 
the system of volume Vq : 



C({nf}) = M-C({n}), 



(1) 



being M a square matrix with coefficients that depend 
only on V and Vq. 

For example, for the Lennard-Jones potential, u(r) = 
4e[((7/V) 12 — (cr/r) 6 ], the vector C can be chosen with 
only two components, Cq and C\ 



Co({n}) = E ( ^ 



<i,3> 



12 



(2a) 



<Wi»=E (£) • 



<M> 



(2b) 



satisfying the two properties mentioned above. 

For an arbitrary potential it may not be possible to 
find variables C n with a given volume scaling as in the 
Lennard-Jones case. However, there is always a method 
based on volume expansions that we describe next. We 
define the coefficient C n from the volume derivatives of 
E{{V/V o y>*{n}): 



Cn({n}) 



d n E((V/V ) 1 / 3 {r i }) 



dV r - 



(3) 



The two properties are fulfilled since the energy of a given 
configuration is E({fi}) — Co({fi}) and the series expan- 
sion, 



c n (m)=f:f^(v-voy- n , 



(4) 



provide the linear relation ([[]). However the vector C 
has an infinite number of components. In practical nu- 
merical work the above expansion needs to be stopped 
at a sufficiently high-order. As it will be seen below the 
approximation introduced can be controlled either by in- 
creasing the order of the approximation or by combining 
simulations at closer densities. 

We denote the density of states with variables C({fif}) 
in some neighborhood of c(V) for a system of volume V, 
by £l(c(V),V). This quantity can be obtained from a 
phase space integration, 

Q(c(n V) = I dr t r... dr N i 6(C({n/}) - c(V)). (5) 
Jv 

Changing the variables of integration, rj = 
{V/VoY^fi and noting relation (|l|) between the vector 
C({fil}) and C({fi}) we can write, 



n(c(v),v) 



N 



V ( , 



dn...df N S(M(C({n})-c(V ))), (6) 



where 



c(V) = M • c(V ). 



(J) 



Using the property of the Dirac delta function, 
S(M(C({r t }) - c(V))) = S(C({r t }) - S(V))/\M\, where 
M| is the determinant of the matrix M, we see that the 
densities of states at different volumes are related by 



V 



N, 



n(c(V),V) dc(V) = (—Y"n(c(V ),V ) dc(V ), (8) 
Vo 

with c(V) and c(Vo) related by equation (0). 



Suppose that we perform several Monte Carlo simula- 
tions at inverse temperature Pi and volume Vi, 1 < i < R. 
Each simulation provides an estimate of the density of 

states: 

n(aw),vo sc(v t ) » ^wwwdffl), (9 ) 

where E(c(Vi)) is the potential energy of the system, 
hi(c(Vi)) is the histogram of the variables C n measured 
in the simulation i, Sc(Vi) is the histogram bin size, Mi is 
the number of measures and /j is the Helmholtz free en- 
ergy at inverse temperature /% and volume Vi . The values 
of fi are not known by now but will be self-consistently 
determined later. 

We use (H) to relate the density of states at volume 
V to the density of states estimated at the simulation 
volume by the above equation, 

n(c(v),v) Sc(v) « f^y e*WW))-f<) h(#yi)) _ 



Mi 



(10) 



The estimates of the density of states given by each of 
the R simulations are now combined |2| , 

Cl(c(V),V) Sc{V) = 

£« (^^«-«, (id 

i=l,R \ */ * 

assigning to each of them a weight pi. The normalized 
(^2i—iPi — 1) weights are obtained from the condition of 
minimization of the statistical uncertainty on the density 
of states, 



<ra(c(F), v) = n 2 {c(v), v) - n(c(v), v) . (12) 

The number of measures in each bin of the histogram is a 
random variable. Neglecting the correlations between the 
measures and using the independence of the simulations 
we have, 



hiWViWMVj)) - hMV)) hMVj)) « hiWViVSij. 

(13) 
The result for the weights is: 

-1 =e 0i(E(S(Vi))-fi) 

xt(^)(vy^ mm ^" 1 (u) 

The partition function at inverse temperature /3 and vol- 
ume V is thus: 



*ow = £E 



hMVi)) e-mzy)) 



^tlEliMdV/Vr e-A^(c(V,))-/«) 

(15) 



and the canonical average of any function f(c(V)) is 

1 



</> = 



Z((3,V) 



k?)hi E^ 1 M I (V I /V)VA(W)H)' 

where X)c(v) ^ s a sum over bins in the multidimensional 
c space. 

For the expansion (Q), the system pressure, P(/3, V), is 
obtained directly from C\: 



p{P,v) = lL-{c 1 ) 



(17) 



It is clear that in the actual calculations there is no need 
to compute the histograms. Denoting by c*> J the measure 
J (1 < J < Mi) in the simulation number i we have: 



(/> = 



1 



Z(J3,V) 

R M, 



f{^{V))e-P E ^ 3( y» 



=w=J Ei=i Mj^/VJ^c-A^C^W))-/!) 



, (18) 



where 



Z(f3,V) = 

R Mi 

EE 



-/3E(c^(V0) 



^!~i EZiMiiVi/V^e-ME^Hv^-f,) 



(19) 



One should remark that in the above expression the 
values <?' J (V;), with / = i are measured while the corre- 
sponding coefficients with I ^ i as well as c* J (V) are com- 
puted from the measured values using equation (m). The 
free energies fa are self-consistently obtained from the 
conditions fa = — (3^ 1 \nZ(/3i, Vi) and by setting fa = 0. 
Thus, we are able to compute free energies relative to 
some thermodynamic state /3i, V\. 



III. APPLICATION TO THE LENNARD-JONES 
FLUID 

We measure values of C every 10 MCS/-/V and the sim- 
ulation lengths were 10 5 MCS/N. The value of the cut-off 
radius was always equal to half the side of the simulation 
box. Standard long range corrections were added to the 
measured values at the end of the simulation. 

We first considered the choice (0) for the vector C. 
Two sets of simulations, with 108 particles, were done 
at two reduced temperatures, Tj* = 1.15 , and T 2 * = 
1.3. For the first temperature we made 40 simulations at 
equally spaced reduced densities: 0.02 < p* < 0.8. For 
the second temperature we made simulations at densities 
p* = 0.02 * LI* -1 , 1 <i <40. 



Every pair of simulations close in density were com- 
bined using the proposed method to obtain results for 
densities in between the two simulations and for a given 
range of temperatures (above and below the simulation 
temperature). From the free energy values obtained we 
built the volume and temperature dependence of the free 
energy. In Fig. \j\ we show the free energy as a func- 
tion of volume per particle at four different temperatures 
T* = 1.0, 1.15, 1.3 and 1.45. In this figure we also com- 
pare the extrapolations obtained from each of the two 
sets of simulations made at T* = 1.15 and T 2 * = 1.3. 
The curves from these two simulations are nearly coinci- 
dent and they are not distinguishable in the figure. For 
T* = 1.15 we also show the common tangent straight 
line at the liquid and gas coexisting phases. The dou- 
ble tangent construction allows us to find the volumes of 
the coexisting phases at each temperature. A new set of 
simulations with 256 particles at T* = 1.3 and densities 
p* = 0.1 * 1.047'- 1 , 1 < i < 40 was also done. 




FIG. 1. Relative Helmholtz free energy as a func- 

tion of volume per particle at four different temperatures 
T* = 1.0, 1.15, 1.3 and 1.45. The results obtained from the 
two simulation sets at temperatures T{ — 1.15 and T 2 * = 1.3 
are plotted and they are not distinguishable. We also show 
the double tangent straight line for T* = 1.15. 

In Fig we show the phase diagram computed from 
each set of simulations. 




FIG. 2. Liquid-Gas phase diagram of the three dimen- 
sional Lennard- Jones model. Solid line: Simulation tempera- 
ture T* = 1.3 and 108 particles; Dashed line: T* = 1.15 and 
108 particles; Dotted line: T* = 1.3 and 256 particles; Circles 
and Triangles are Gibbs Ensemble results from reference lid] 
with 500 particles and 300 particles respectively. 



In order to ascertain the usefulness of volume expan- 
sion method (0) we also made simulations where the vari- 
ables, C n , < n < 5, were measured. This set of 
simulations was done for a system of 108 particles at 
T* = 1.3 and at the same densities chosen before. In Fig 

we show the convergence of the results for the relative 
free energy as a function of volume for a temperature 
T* = 1.15 as we include an increasing number of coeffi- 
cients in the expansion (|4|). The curves obtained with 4, 
5 and 6 coefficients are nearly coincident and agree with 
results obtained from the choice based on equation (g). 
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FIG. 3. Convergence of the free energy with an increasing 
number of coefficients in the expansion based in Eq. ffl). 
Results with 1 to 6 coefficients correspond respectively to 
long-dashed, dot dashed, dot-dot-dashed, short-dashed, dot- 
ted and solid curves. The curves obtained with 4, 5 and 6 
coefficients are nearly coincident. 



IV. CONCLUSIONS 

We have proposed a method which allows simultane- 
ous extrapolations in volume and temperature based on 
the multiple histogram method. An arbitrary number of 
Monte Carlo simulations made in the canonical ensemble 
can be combined providing improved estimates of ther- 
modynamic properties. We show test bed results on the 
three-dimensional Lennard- Jones system which confirm 
that the method works well and that the volume expan- 
sion scheme based on equation (^) can be used with a 
good control of the approximations involved. Calcula- 
tion of relative Helmholtz free energy coupled with the 
double tangent construction allows an efficient determi- 
nation of the phase diagram. The method is general and 
can be applied to interaction potentials that do not have 
a simple scaling with system volume. 



ACKNOWLEDGEMENTS 

We thank Professors E. J. S. Lage, S. K. Mendi- 
ratta, T. P. Gasche and J. M. Pacheco for a criti- 
cal reading of the manuscript. A. L. Ferreira grate- 
fully acknowledges J. M. Pacheco and J. P. Ra- 
malho for fruitful discussions. This work was partially 
supported by the projects PRAXIS/2/2.1/299/94 and 
PRAXIS/2/2. 2/FIS/302/94. 



[1] A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 

61, 2635 (1988). 
[2] A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 

63, 1195 (1989). 
[3] R. H. Swend sen, Physica A 194 , 53-62 (1993). 
[4] E. Marinari, |cond-mat/961201C| (1996). 
[5] G. M. Torrie and J. P. Valleau, J. Chem. Phys. 66, 1402 

(1977). 
[6] J. P. Valleau, J. Comput. Phys. 96, 193 (1991). 
[7] N. V. Brilliantov and J. P. Valleau, J. Comput. Phys. 

108, 1115 (1998). 
[8] B. A. Berg and T. Neuhaus, Phys. Rev. Lett. 68, 9 

(1992). 
[9] W. Janke, in Computer Simulations in Condensed Matter 
Physics VII, edited by D. P. Landau, K. K. Mon and H. 
B. Schuttler (Springer-Verlag, Heidelberg, 1994). 
[10] A. Z. Panagiotopoulos, Mol. Phys. 61, 813-826 (1987). 
[11] A. Z. Panagiotopoulos, Mol. Simulation. 9, 1 (1992). 
[12] J. M. Caillol, J. Chem. Phys 109, 4885 (1998). 
[13] N. B. Wilding, Phys. Rev. E 52, 602 (1995). 
[14] A. D. Bruce and N. B. Wilding, Phys. Rev. Lett. 68, 193 

(1992). 
[15] A. P. Lyubartsev et. al., J. Chem. Phys. 96, 1776 (1992). 



[16] F. A. Escobedo and J. de Pablo, Mol. Phys. 87, 347 

(1996). 
[17] E. J. Meijer et. al., J. Chem. Phys. 92, 7570 (1990). 
[18] A. L. Ferreira, J. M. Pacheco and J. P. Ramalho, work 

in progress. 
[19] D. A. Kofke, J. Chem. Phys. 98, 4149 (1993). 
[20] M. Mehta and D. A. Kofke, Mol. Phys. 86, 139 (1995). 
[21] F. A. Escobedo, J. Chem. Phys. 108, 8761 (1998). 



